#---------------------------------------------------------------------------------#
### information: 
#---------------------------------------------------------------------------------#
### authors: björn bremer (bremer@mpifg.de), reto bürgisser (buergisser@ipz.uzh.ch)
### last updated: 22 October 2021

#---------------------------------------------------------------------------------#
### description
#---------------------------------------------------------------------------------#

### this file replicates the analysis of the following paper:
### "Public Opinion on Welfare State Recalibration in Times of Austerity"
### this file focuses on the results from the SPLIT SAMPLE EXPERIMENT shown in the MAIN TEXT

#---------------------------------------------------------------------------------#

set.seed(1234567)

#--------------------------------------------------------
### load libraries ###
#--------------------------------------------------------

library(grid)
library(gridExtra)
library(ggplot2)
library(ggpubr)
library(margins)
library(plyr)
library(dplyr)

#--------------------------------------------------------
### load data ###
#--------------------------------------------------------

load("df_split_psrm_replication.Rdata") # load the data file

#--------------------------------------------------------
### figure 1: support for social investment and consumption ###
#--------------------------------------------------------

## social consumption
bar_sc <- ggplot(df_split, aes(x=sc)) +
      geom_bar(aes(group=factor(0),y = ..prop..), width = 0.9) + ylim(0,0.41) + 
      theme_bw() + 
      scale_x_continuous(breaks=seq(0,10,1)) +
      ylab("Proportion") +
      xlab("Support for social consumption (0-10)") +
      ggtitle("Social Consumption") + theme(plot.title = element_text(hjust = 0.5, face = "bold")) 
bar_sc

## social investment
bar_si <- ggplot(df_split, aes(x=si)) +
      geom_bar(aes(group=factor(0),y = ..prop..), width = 0.9) + ylim(0,0.41) + 
      theme_bw() + 
      scale_x_continuous(breaks=seq(0,10,1)) +
      ylab("Proportion") +
      xlab("Support for social investment (0-10)") +
      ggtitle("Social Investment") + theme(plot.title = element_text(hjust = 0.5, face = "bold")) 
bar_si

## combine the graphs and save the resulting figure
bar_motivation <- grid.arrange(bar_sc,bar_si, ncol = 2)
ggsave(bar_motivation , width = 20, height = 12, units = c("cm"), file ="figure1.eps")

#--------------------------------------------------------
### figure 2: average support for spending increases by treatment and model, pooled ###
#--------------------------------------------------------

### calculate ols regressions

## education
edu_1 <- lm(edu ~ as.factor(split), data = df_split)
edu_2 <- lm(edu ~ as.factor(split) + age + female + children_u10 + education_c + occup_gh + workctr + party_fam_rec, data = df_split)
edu_3 <- lm(edu ~ as.factor(split) + age + female + children_u10 + education_c + occup_gh + workctr + party_fam_rec + country, data = df_split)

## almp
almp_1 <- lm(almp ~ as.factor(split), data = df_split)
almp_2 <- lm(almp ~ as.factor(split) + age + female + children_u10 + education_c + occup_gh + workctr + party_fam_rec, data = df_split)
almp_3 <- lm(almp ~ as.factor(split) + age + female + children_u10 + education_c + occup_gh + workctr + party_fam_rec + country, data = df_split)

## childcare
childcare_1 <- lm(childcare ~ as.factor(split), data = df_split)
childcare_2 <- lm(childcare ~ as.factor(split) + age + female + children_u10 + education_c + occup_gh + workctr + party_fam_rec, data = df_split)
childcare_3 <- lm(childcare ~ as.factor(split) + age + female + children_u10 + education_c + occup_gh + workctr + party_fam_rec + country, data = df_split)

### predict values

## education
pre_edu_1 <- prediction(edu_1, df_split)
pre_edu_2 <- prediction(edu_2, df_split)
pre_edu_3 <- prediction(edu_3, df_split)

# calculate average predicted values by split
sum_edu_1 <- ddply(pre_edu_1, c("split"), summarise,
                   mean = mean(fitted, na.rm=TRUE),
                   mean_se = mean(se.fitted, na.rm=TRUE))
sum_edu_1$model <- "M1"

sum_edu_2 <- ddply(pre_edu_2, c("split"), summarise,
                   mean = mean(fitted, na.rm=TRUE),
                   mean_se = mean(se.fitted, na.rm=TRUE))
sum_edu_2$model <- "M2"

sum_edu_3 <- ddply(pre_edu_3, c("split"), summarise,
                   mean = mean(fitted, na.rm=TRUE),
                   mean_se = mean(se.fitted, na.rm=TRUE))
sum_edu_3$model <- "M3"

# combine all estimates in a single df
sum_edu_all <- rbind(sum_edu_1, sum_edu_2, sum_edu_3)

## ALMP
pre_almp_1 <- prediction(almp_1, df_split)
pre_almp_2 <- prediction(almp_2, df_split)
pre_almp_3 <- prediction(almp_3, df_split)

# calculate average predicted values by split
sum_almp_1 <- ddply(pre_almp_1, c("split"), summarise,
                    mean = mean(fitted, na.rm=TRUE),
                    mean_se = mean(se.fitted, na.rm=TRUE))
sum_almp_1$model <- "M1"

sum_almp_2 <- ddply(pre_almp_2, c("split"), summarise,
                    mean = mean(fitted, na.rm=TRUE),
                    mean_se = mean(se.fitted, na.rm=TRUE))
sum_almp_2$model <- "M2"

sum_almp_3 <- ddply(pre_almp_3, c("split"), summarise,
                    mean = mean(fitted, na.rm=TRUE),
                    mean_se = mean(se.fitted, na.rm=TRUE))
sum_almp_3$model <- "M3"

# combine all estimates in a single df
sum_almp_all <- rbind(sum_almp_1, sum_almp_2, sum_almp_3)

### create the graphs
## childcare
pre_childcare_1 <- prediction(childcare_1, df_split)
pre_childcare_2 <- prediction(childcare_2, df_split)
pre_childcare_3 <- prediction(childcare_3, df_split)

# calculate average predicted values by split
sum_childcare_1 <- ddply(pre_childcare_1, c("split"), summarise,
                         mean = mean(fitted, na.rm=TRUE),
                         mean_se = mean(se.fitted, na.rm=TRUE))
sum_childcare_1$model <- "M1"

sum_childcare_2 <- ddply(pre_childcare_2, c("split"), summarise,
                         mean = mean(fitted, na.rm=TRUE),
                         mean_se = mean(se.fitted, na.rm=TRUE))
sum_childcare_2$model <- "M2"

sum_childcare_3 <- ddply(pre_childcare_3, c("split"), summarise,
                         mean = mean(fitted, na.rm=TRUE),
                         mean_se = mean(se.fitted, na.rm=TRUE))
sum_childcare_3$model <- "M3"

# combine all estimates in a single df
sum_childcare_all <- rbind(sum_childcare_1, sum_childcare_2, sum_childcare_3)

## education
sum_edu_all$split <- factor(sum_edu_all$split, levels = c(1,3,2,4)) # re-order the graph
plot_pre_edu <- ggplot(sum_edu_all,aes(as.factor(split), mean, colour = model)) +
      geom_errorbar(aes(ymin=mean-(1.96*mean_se), ymax=mean+(1.96*mean_se)), width=0, position=position_dodge(width=0.5)) + 
      geom_point(aes(y=mean, shape = model), position = position_dodge(width=0.5), size=2) +
      scale_y_continuous(limits=c(3, 9), breaks = round(seq(3, 9, by = 1), 1)) +
      ylab("Mean support for spending") +
      theme_bw() +
      theme(axis.text.x = element_text(angle = 45, hjust = 1)) +  
      scale_x_discrete(name = "Treatment", breaks = c(1,2,3,4), c("Control", "Lower ALMP", "Lower pension", "Lower childcare")) +      
      scale_shape_discrete(name = "Models", labels = c("No covariates", "Covariates", "Covariates + country-FE")) +
      scale_colour_grey(name = "Models", labels = c("No covariates", "Covariates", "Covariates + country-FE"), start = 0.2, end = 0.8, aesthetics = "colour") +
      ggtitle("Education") +theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
      coord_fixed(ratio = 0.8)
plot_pre_edu

## active labor market policy
plot_pre_almp <- ggplot(sum_almp_all,aes(as.factor(split), mean, colour = model)) +
      geom_errorbar(aes(ymin=mean-(1.96*mean_se), ymax=mean+(1.96*mean_se)), width=0, position=position_dodge(width=0.5)) + 
      geom_point(aes(y=mean, shape = model), position = position_dodge(width=0.5), size=2) +
      scale_y_continuous(limits=c(3, 9), breaks = round(seq(3, 9, by = 1), 1)) +
      ylab("Mean support for spending") +
      theme_bw() +
      theme(axis.text.x = element_text(angle = 45, hjust = 1)) +  
      scale_x_discrete(name = "Treatment", breaks = c(1,2,3,4), c("Control", "Lower pension", "Lower childcare", "Lower PLMP")) +      
      scale_shape_discrete(name = "Models", labels = c("No covariates", "Covariates", "Covariates + country-FE")) +
      scale_colour_grey(name = "Models", labels = c("No covariates", "Covariates", "Covariates + country-FE"), start = 0.2, end = 0.8, aesthetics = "colour") +
      ggtitle("Active Labor Market Policy") +theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
      coord_fixed(ratio = 0.8) 
plot_pre_almp

## childcare
plot_pre_childcare <- ggplot(sum_childcare_all,aes(as.factor(split), mean, colour = model)) +
      geom_errorbar(aes(ymin=mean-(1.96*mean_se), ymax=mean+(1.96*mean_se)), width=0, position=position_dodge(width=0.5)) + 
      geom_point(aes(y=mean, shape = model), position = position_dodge(width=0.5), size=2) +
      scale_y_continuous(limits=c(3, 9), breaks = round(seq(3, 9, by = 1), 1)) +
      ylab("Mean support for spending") +
      theme_bw() +
      theme(axis.text.x = element_text(angle = 45, hjust = 1)) +  
      scale_x_discrete(name = "Treatment", breaks = c(1,2,3,4), c("Control", "Lower child bf.", "Lower pension", "Lower PLMP")) +      
      scale_shape_discrete(name = "Models", labels = c("No covariates", "Covariates", "Covariates + country-FE")) +
      scale_colour_grey(name = "Models", labels = c("No covariates", "Covariates", "Covariates + country-FE"), start = 0.2, end = 0.8, aesthetics = "colour") +
      ggtitle("Childcare") +theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
      coord_fixed(ratio = 0.8) 
plot_pre_childcare

## combine the graphs and save the resulting figure 
socinv_split_1row <- ggarrange(plot_pre_edu, plot_pre_almp, plot_pre_childcare, ncol=3, nrow=1, common.legend = TRUE, legend="bottom")
ggsave(socinv_split_1row, width = 20, height = 12, units = c("cm"), file ="figure2.eps")

#--------------------------------------------------------
### figure 3: share of supporters for spending increases by treatment ###
#--------------------------------------------------------

## education

# create new education binary support variable 
df_split_bin <- df_split[which (!is.na(df_split$split)),]

df_split_bin$edu_bin <- ifelse(df_split$edu <= 5, 0,
                               ifelse(df_split$edu > 5, 1, NA))

df_split_bin$edu_bin <- df_split_bin$edu_bin*100

# calculate means and se overall
mean_edu_bin <- ddply(df_split_bin, c("split"), summarise,
                      N = sum(!is.na(edu_bin)),
                      mean = mean(edu_bin,  na.rm=TRUE),
                      sd = sd(edu_bin,  na.rm=TRUE),
                      se = sd / sqrt(N)
)

mean_edu_bin$split <- as.factor(mean_edu_bin$split)

# create the graph
barplot_mean_edu_bin <- ggplot(mean_edu_bin, aes(x = split, y = mean)) +
      geom_hline(yintercept = 50, colour = "gray50") + 
      geom_bar(stat="identity", width=0.6) +
      geom_errorbar(aes(ymin=mean-(1.96*se), ymax=mean+(1.96*se)), width=0.3, color="grey20") + 
      scale_y_continuous(limits=c(0,100), breaks = round(seq(0, 100, by = 10), 1)) +
      ylab("Support for education spending (in %)") +
      scale_x_discrete(name = "Treatment", limits = c("1", "3", "2", "4"), labels = c("Control", "Lower pension", "Lower ALMP","Lower childcare")) +  
      theme_bw() +
      theme(axis.text.x = element_text(angle = 45, hjust = 1))  +
      ggtitle("Education") +theme(plot.title = element_text(hjust = 0.5, face = "bold")) 
barplot_mean_edu_bin

## almp

# create new almp binary support variable 
df_split_bin <- df_split[which (!is.na(df_split$split)),]

df_split_bin$almp_bin <- ifelse(df_split$almp <= 5, 0,
                                ifelse(df_split$almp > 5, 1, NA))

df_split_bin$almp_bin <- df_split_bin$almp_bin*100

# calculate means and se overall
mean_almp_bin <- ddply(df_split_bin, c("split"), summarise,
                       N = sum(!is.na(almp_bin)),
                       mean = mean(almp_bin,  na.rm=TRUE),
                       sd = sd(almp_bin,  na.rm=TRUE),
                       se = sd / sqrt(N)
)

mean_almp_bin$split <- as.factor(mean_almp_bin$split)

# create the graph
barplot_mean_almp_bin <- ggplot(mean_almp_bin, aes(x = split, y = mean)) +
      geom_hline(yintercept = 50, colour = "gray50") + 
      geom_bar(stat="identity", width=0.6) +
      geom_errorbar(aes(ymin=mean-(1.96*se), ymax=mean+(1.96*se)), width=0.3, color="grey20") + 
      scale_y_continuous(limits=c(0,100), breaks = round(seq(0, 100, by = 10), 1)) +
      ylab("Support for ALMP spending (in %)") +
      scale_x_discrete(name = "Treatment", limits = c("1", "4", "2", "3"), labels = c("Control", "Lower PLMP", "Lower pension", "Lower childcare")) +  
      theme_bw() +
      theme(axis.text.x = element_text(angle = 45, hjust = 1))  +
      ggtitle("Active Labor Market Policy") +theme(plot.title = element_text(hjust = 0.5, face = "bold")) 
barplot_mean_almp_bin

## childcare

# create new education binary support variable 
df_split_bin <- df_split[which (!is.na(df_split$split)),]

df_split_bin$childcare_bin <- ifelse(df_split$childcare <= 5, 0,
                                     ifelse(df_split$childcare > 5, 1, NA))

df_split_bin$childcare_bin <- df_split_bin$childcare_bin*100 

# calculate means and se overall
mean_childcare_bin <- ddply(df_split_bin, c("split"), summarise,
                            N = sum(!is.na(childcare_bin)),
                            mean = mean(childcare_bin,  na.rm=TRUE),
                            sd = sd(childcare_bin,  na.rm=TRUE),
                            se = sd / sqrt(N)
)

mean_childcare_bin$split <- as.factor(mean_childcare_bin$split)

# creat the graph
barplot_mean_childcare_bin <- ggplot(mean_childcare_bin, aes(x = split, y = mean)) +
      geom_hline(yintercept = 50, colour = "gray50") + 
      geom_bar(stat="identity", width=0.6) +
      geom_errorbar(aes(ymin=mean-(1.96*se), ymax=mean+(1.96*se)), width=0.3, color="grey20") + 
      scale_y_continuous(limits=c(0,100), breaks = round(seq(0, 100, by = 10), 1)) +
      ylab("Support for childcare spending (in %)") +
      scale_x_discrete(name = "Treatment", labels = c("Control", "Lower child bf.", "Lower pension", "Lower PLMP")) +  
      theme_bw() +
      theme(axis.text.x = element_text(angle = 45, hjust = 1))  + 
       ggtitle("Childcare") +theme(plot.title = element_text(hjust = 0.5, face = "bold")) 
barplot_mean_childcare_bin

## combine the graphs and save the resulting figure 
sociinv_split_bar <- ggarrange(barplot_mean_edu_bin, barplot_mean_almp_bin, barplot_mean_childcare_bin, ncol = 3, nrow = 1)
ggsave(sociinv_split_bar, width = 25, height = 11, units = c("cm"), file ="figure3.eps")
